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Abstract 

The size and geometry of the prostate are known to be pivotal quanti- 
ties used by clinicians to assess the condition of the gland during prostate 
cancer screening. As an alternative to palpation, an increasing number 
of methods for estimation of the above-mentioned quantities are based on 
using imagery data of prostate. The necessity to process large volumes of 
such data creates a need for automatic segmentation tools which would 
allow the estimation to be carried out with maximum accuracy and effi- 
ciency. In particular, the use of transrectal ultrasound (TRUS) imaging in 
prostate cancer screening seems to be becoming a standard clinical prac- 
tice due to the high benefit-to-cost ratio of this imaging modality. Unfor- 
tunately, the segmentation of TRUS images is still hampered by relatively 
low contrast and reduced SNR of the images, thereby requiring the seg- 
mentation algorithms to incorporate prior knowledge about the geometry 
of the gland. In this paper, a novel approach to the problem of segment- 
ing the TRUS images is described. The proposed approach is based on 
the concept of distribution tracking, which provides a unified framework 
for modeling and fusing image-related and morphological features of the 
prostate. Moreover, the same framework allows the segmentation to be 
regularized via using a new type of "weak" shape priors, which minimally 
bias the estimation procedure, while rendering the latter stable and ro- 
bust. The value of the proposed methodology is demonstrated in a series 
of both in silico and in vivo experiments. 
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1 Introduction 



Ultrasound imaging has long become an integral part of modern health care 
due to its superiority over many alternative imaging modalities in terms of its 
cost-to-benefit ratio. In particular, the properties of ultrasound imaging be- 
ing practically harmless, mobile, readily accessible, and cost efficient has made 
it the modality of choice in many clinical settings. Unfortunately, the above 
advantages have always been counterbalanced by the relatively low quality of 
ultrasound images as compared to X-ray CT or MRI scans. As a result, there 
has always been a need for image processing tools which, when applied to ultra- 
sound images, could produce valuable outcomes under the conditions of reduced 
resolution and contrast [1]. 

According to recent global cancer statistics, prostate cancer constitutes the 
fifth most common type of cancer in the world and the second most common 
in men [2]. The medical treatment cost of prostate cancer runs into millions 
of dollars annually, while the emotional cost to the patients and the members 
of their families is incalculable. This is why early diagnosis of prostate cancer 
through systematic screening has become the main tool of prevention of the dis- 
ease. Even though the current practice of screening is still based on palpation, it 
is believed that more reliable diagnostic outcome is possible via using minimally 
invasive procedures, such as transrectal ultrasound imaging (TRUS) [3]. More- 
over, since the size and shape of the prostate are among the main quantities 
used by clinicians to identify the presence of abnormal developments within the 
gland [4], automatic segmentation of ultrasound prostate images is currently 
recognized as a pivotal component of the TRUS-based diagnosis of prostate 
cancer [5-9]. 

In the case of TRUS imaging, the existing methods of prostate segmenta- 
tion encompass a number of different strategies. Some basic solutions to the 
problem are based on utilizing edge-related features pertaining to the prostate 
boundary [5,7, 10]. For example, the method of [10] prompts the user to man- 
ually input a set of initial points on the boundary of the gland. These points 
are used to initialize a discrete dynamic contour [11], which is subsequently 
driven towards strong local edges. Alternatively, the approach of [5] detects the 
prostate boundary using the Canny edge detector applied to the TRUS images 
enhanced by means of anisotropic diffusion filtering. The resulting edge maps 
are then superimposed over the original images as a visual guide to subsequent 
manual delineation. Note that, even though this method has not been devised 
as a fully automatic procedure, it could still be useful to reduce both intra- and 
inter-observer variances of prostate delineation. To improve the robustness of 
prostate segmentation in the case of poorly observable edges, a two step segmen- 
tation procedure was proposed in [7]. At the first stage of this method, a coarse 
segmentation is obtained based on texture-related features of the TRUS images, 
followed by forcing the segmentation boundary to converge to local image edges 
during the second stage. Unfortunately, this method (as well as all the methods 
mentioned before) has a drawback of being prone to the errors caused by the 
existence of spurious (noise-induced) edges and shadowing artifacts. 
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Another classical approach to the problem of prostate segmentation takes 
advantage of the tools of supervised [12,13] and unsupervised [14] machine learn- 
ing. Specifically, the method of [12] employs a feed-forward neural network to 
carry out the segmentation, which attempts to capture the visual cues used by 
experts to classify image pixels as either inside or outside of the prostate region. 
In the recent approach of [13], the image features used by a support vector 
machine classifier are defined via projecting the TRUS images onto a basis of 
multi- wavelets. On the side of unsupervised learning, a probabilistic clustering 
procedure detailed in [14] segments the prostate region based on the features 
obtained using a specially designed filter bank [15]. Despite the strong theoreti- 
cal foundation of the classification-based segmentation algorithms, they remain 
prone to over-segmentation because of their disregard for the morphological 
features of the prostate gland. 

Recent efforts on prostate segmentation have been extended to use the 
information provided by shape priors to regularize the process of segmenta- 
tion [6,8,16,17]. Thus, for example, the method in [6] models the prostate 
boundary as a deformable super-ellipse. Subsequently, a maximum-a-posteriori 
estimation framework is used to find a contour (i.e. a boundary of the prostate) 
that closely matches the prior shape model on one hand, and coincides with 
strong image gradients on the other. The same super-ellipse model is employed 
in [16] where texture-related features are used instead of edge-related ones. An 
obvious limitation of the above approaches stems from the restrictive nature 
of the prior model used, which can only approximate symmetric shapes. More 
general assumptions regarding the shape of prostate are used in [17], which uses 
a Gabor filter bank to describe the prostate boundary. Subsequently, a hier- 
archical shape deformation procedure moves an initial segmentation towards 
the location of the true boundaries of the gland. The implementation of this 
method, however, is hindered by a number of initialization procedures such as a 
normalization of the shape of the ultrasound transducer in use, and a geometric 
transformation of the TRUS images to align them with a prior shape model. 
Moreover, the use of a greedy energy-minimization algorithm to search for an 
optimal prostate boundary leads to relatively large computation times, which 
makes the above method impractical to apply to large data sets. 

Among other methods for prostate segmentation, one can mention the semi- 
automatic procedure of [18] which is based on the method of graph cuts, and 
the approach in [19] which takes advantage of a Kalman state estimator [20] to 
track the prostate boundary starting from a seed location. 

In the current paper, a different approach to the problem of segmentation of 
TRUS images is introduced. As opposed to the previous algorithms described 
above, the proposed method is exceptional for it concurrently fulfills a number 
of essential objectives, viz. 

1. Generality: The algorithm fuses information from a range of different 
sources, viz. edge-related, texture-related, and shape-related features are 
all used to estimate an optimal boundary of the prostate gland. 

2. Adaptability: The method allows using any number of arbitrary texture- 
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related features of the TRUS images such as their intensity, local polyno- 
mial moments, multiresolution coefficients, local entropies, etc. Moreover, 
the method can also accommodate the use of any number of shape de- 
scriptors (subject to some transformation-invariance requirements to be 
discussed below). 

3. Consistency: The method takes advantage of the concept of distribution 
tracking [21], which allows applying the same methodology to handle both 
photometric and morphological features of the prostate, thereby providing 
a unified framework for the segmentation problem at hand. 

Finally, it is important to note that integrating shape prior models into 
probabilistic image segmentation has never been a trivial task, with a spectrum 
of different approaches proposed to this end [22]. In most of the cases, the an- 
tecedent knowledge on plausible configurations of an optimal shape is "encoded" 
in the form of prior probabilities which can describe the distribution of, e.g., the 
control points [23] or level-set functions [24] related to the shape representation. 
Such methods have proven particularly useful in the cases when the objects to 
be segmented appear to be partially observable due to, e.g., occlusions. Under 
these circumstances, the use of elaborate shape priors seems to be well justified, 
even when a sizable bias is introduced into the estimated shapes. For the case 
at hand, however, using such priors would not be a proper choice, since the 
prostates are rarely (if ever) occluded in TRUS imaging. In this situation, the 
prior information should be restrictive just enough to regularize the process of 
segmentation, which should still be dominated by the information contained in 
the acquired images. 

To achieve the above goal, a new type of "weak" shape priors is introduced 
in this paper. The prior information is encapsulated in the form of a probability 
density function (pdf) of an intrinsic geometric parameter (or a set thereof) of 
the true prostate boundary. Subsequently, the segmentation boundary is forced 
to evolve into a configuration whose empirical pdf of the same geometric param- 
eter^) closely matches the model pdf as assessed by the Bhattacharyya distance 
measure. In this sense, the present paper extends the approaches of [21,25] to 
tracking both morphological and photometric features of the object of interest. 
Moreover, it is proven experimentally that the proposed methodology can be 
effectively used for robust and accurate segmentation of TRUS images under the 
conditions of low image contrast and poorly observable prostate boundaries. 

The rest of the paper is organized as follows. In Section II, a method for 
prostate image segmentation via tracking of both texture-related and morpho- 
logical features is detailed. A number of essential technical considerations are 
provided in Section III, while Section IV presents a series of both in silico and 
in vivo experiments. Finally, Section V summarizes the paper with a discussion, 
conclusions, and an outline of our future research. 
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2 Tracking of distributions 



2.1 Tracking texture-related features 

Let u(x) be an ultrasound image which we consider to be a scalar- valued function 
defined over a subset ft of R 2 . To render our discussion general, let M be a map 
which transforms it(x) into its corresponding feature image /(x) = .M[it(x)], 
which can be viewed as a vector-valued image (i.e. /(x) : f2 — > R d ) with 
its components being the features pertaining to it(x). Thus, for each x e il, 
J(x) may consist of the values of the original image u(x), its multiresolution 
version [26], some phase-related and similarity measures [27], etc. 

Let il t (with the subscript t standing for "target" ) be a subset of over which 
the object of interest (i.e. prostate) is supported. The segmentation algorithm 
proposed in this paper is based on the following two assumptions. First, for 
each x 6 f2 t , the d components of /(x) are assumed to be independent random 
variable^] Second, it is assumed that one is provided with reliable estimates 
of the probability density functions (pdf) of the image features pertaining to 
the prostate region £l t . Note that, even though the second assumption could 
seem impractical, it is actually easy to attain due to the availability of vast 
depositories of manually segmented prostate images. These images can readily 
be used to extract the values of the features of interest within the boundaries 
of manually delineated prostates. These values in turn can be used to estimate 
their corresponding pdf's by means of any standard procedure [28]. 

Let -Pt(z), where z £ R d , be the joint pdf of the image features associated 
with the prostate region fl t . Due to the assumption of statistical independence 
of the image features, P(z) can be factorized as P(z) = nf=i Pt( z k), with 
Pt(zk) being the pdf of the fc-th image feature z k € R. Now, let C be an 
arbitrary subset of f2 (which is not equal to fl t , in general). Then, for a given 
/(x) one can define the empirical probability density of z k observed over Q in as 

f in ^ In^ K ( z k^h^))d X 

p{z k n in ) = — , (i) 

where /^(x) is the fc-th component of /(x), and K is a positive- valued kernel 
function satisfying J K(x)dx — 1. Note that the expression in is nothing else 
but a kernel-based, nonparametric estimate of the pdf of z k [28, 29]. By virtue 
of the statistical independence of z k , the joint (empirical) pdf of z observed over 
flin is given by P(z | Qi n ) = Yit=iP( z k I ^m)- Consequently, one can quantify 
the degree of similarity between the target P*(z) and the empirical P(z | f2j n ) 
densities by means of the Bhattacharyya coefficient B as given by [30-32] 

B(n in )= [ y/p t (x)p(z | n m )dz. (2) 

1 Note that the assumption of independence has been mainly introduced to render the algo- 
rithm computationally feasible, and thus a mitigation of this assumption should be addressed 
in future research. 
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Alternatively, using the assumption of mutual independence of z k , one can 
rewrite ^ as 

d 

B(n in ) = Y[B k (Q in ), (3) 

k=l 

where 

B k (n in) = l^ Pt ( Zk ) p{Zk \n in ) dZk . (4) 

Viewed as a function of flj n , B(f2j n ) takes its values in the interval [0, 1], 
with the maximum attained when pt(z k ) = p(z k | f2i„),Vfc = 1,2, ...,d. In 
particular, B(fli n ) is maximized when the regions £li n and fit coincide with each 
other (meaning f2j n = Q t ). Consequently, the image w(x) can be segmented via 
approximating f2 t by fl in which locally maximizes B(fl in ) [21]. 



2.2 Level-set segmentation paradigm 

It goes without saying that maximizing B(£li n ) as a function of f2, n would result 
in a combinatorial optimization of little practical interest. To find a useful 
segmentation in a consistent manner, a different formulation of the problem is 
needed. In particular, in this paper, fli n is defined implicitly as a subset of Q, 
over which a level-set function y(x) : f2 — > K attains non-positive values, wz. 

fiin = {xenj p(x) < 0} . (5) 

Consequently, using the standard definition of the Heaviside function as Ti (x) = 
(#)+, the Bhattacharyya coefficient in ^ can be redefined as a function of y(x) 
as given by 

S(V>(x)) = [] B»((p(x)) = 1] / VftRK^Mx))^, (6) 

where 

, , , „ / n IT(z fc -/ fe (x))ft(-y(x))rfx 

v?(x)) = r — -— . 7 

Consequently, the problem of segmentation of u(x) amounts to finding an opti- 
mal </?(x) which maximizes i?((/?(x)). 

A (local) maximizer of £?(<£>(x)) can be computed via the standard steepest 
ascent procedure which has the form of a gradient flow defined as 

/__ \ a d<p(x,T) 6B(<p) 
^(x,r) = ^— = — , (8) 

with r being an artificial time parameter (so that </?(x, 0) = v 5 o( x ) is an initial- 
ization of the level-set function), and SB (if) /Sip standing for the first variation 
of B(ip(x)) computed with respect to y(x). Appendix A provides necessary 
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technical details on the derivation of 5B(ip)/8tp, and shows that ^ can be 
rewritten as given by 

<^(x,t) = %(x,t))Vb(x,t), (9) 

where S(-) stands for the Dirac delta function, and S(ip(-)) Vb(-) ■ — * K is a 
velocity field that governs the evolution of the level-set y(x). (In the develop- 
ments below, we suppress the explicit dependency of both level-set functions 
and velocity fields on r, until it is needed to help clarify the derivations). 

It should be noted that the form of the gradient flow in ^ implies that 
the latter depends on the observed data /(x) alone, while disregarding some 
plausible properties of the optimal solution. As a result, maximizing (|6| could be 
too sensitive to measurement noises and/or errors in the data. This sensitivity, 
however, can be alleviated via regularizing the maximization of B(ip(x)) using 
the framework of geodesic active contours [33,34]. In this case, the optimal </?(x) 
is found as 

p*(x) = argsup faB(p(x)) - / <?(x) ||W%(x))|| dx) , (10) 
ip(x) I Jn ) 

where V denotes the operator of gradient, || ■ || is the Euclidean norm, a > is 
a regularization constant]^ and g : — > 1R+ is an edge-detector function (to be 
defined below). As a result, the gradient flow associated with the optimization 



problem in (10) can be shown to be equal to 



p r (x) = %(x)) (a Vs(x) +div (ff( x ) ||v^|x)|| )) ' (U) 

It is worth pointing out that the gradient flow in ( |11[ ) is designed to converge 
to an optimal level-set function, whose related £li„ has its boundary aligned with 
the strong edges of it(x), while the empirical distributions of the image features 
observed over Q in are "aligned" with the model distributions pt{z^) in the sense 
of maximizing B(ip(x)). 



Finally, it should be noted that the gradient flow (111 would change the 
level-set function y(x) over a set of zero measure, if the formal delta function 
S(-) was used in computations. To overcome this "technical" difficulty, it is 
common to extend the numerical support of the level-set evolution via replacing 
S(-) by its smoother version 5 e (-) that could be defined as given by, e.g., [35]: 

{<w= /4(i + «ro). ms. (12) 

I 0, otherwise 

Note that the support of S e is finite and controlled though the user-defined 
parameter e. Throughout the experimental study of this paper, e is set to be 
equal to 2. 



2 In general, a can be a user defined parameter. In the current paper, its value has been 
set to be equal to 0.5. 
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Original prostate image Edge -detector g[x) 




Figure 1: Segmentation of prostate image without shape priors: (A) Original 
image; (B) Edge-detector function g(x); (C) Initial segmentation; (D) Final 
segmentation. 

2.3 A Motivating Example 

In practical computations, the boundary Ti n of 0j„ is used akin to a decision 
boundary that separates the object of interest from its surrounding. This bound- 
ary - conventionally referred to as an active contour - is defined by the zero 
level set of <p( x ), i-e- r in = {x € f2 | y(x) = 0}. In the present case, the active 
contour belongs to the family of geodesic active contours, which are particularly 
useful in the situations when objects of interest have well-defined edges. Unfor- 
tunately, in ultrasound imaging, edges can hardly be considered as reliable fea- 
tures. Due to shadowing artifacts, diffraction, as well as the destructive nature 
of speckle noise, the ultrasound representation of many anatomical structures 
(including prostates) happens to lack continuously defined boundaries. Conse- 
quently, using the segmentation tools dependent on an edge-detection procedure 
is prone to converge to erroneous results in the case of ultrasound imaging. 
The above problem is exemplified in Fig. [I] which shows an original prostate 



8 



image (Subplot A), its corresponding edge-detector function g(x) (Subplot B), 
an initial segmentation (Subplot C), and a final segmentation (Subplot D) as 
obtained by ( |TT] ). Note that in this experiment, the edge-detector function g(x) 
was defined aqfj 

5(x) = i+Aiiw(x)r (13) 

where u(x) stands for a de-speckled version of it(x) computed by the SRAD 
filter [36], and A is a positive constant (set to be equal to 3 in this example). 
Moreover, the Bhattacharyya functional in ^ was defined using ii(x) = u(x) 
and l2(x) = w(x). One can see that the incompleteness of the prostate boundary 
along with the shadowing artifacts have caused the active contour to "leak out" 
at the sides of the prostate region. To alleviate this deficiency, additional infor- 
mation on the expected geometry of prostate boundary should be incorporated 
into the segmentation procedure, as explained next. 

2.4 Tracking of Shape Features using Weak Priors 

To overcome the limitation presented in the example of Fig. [I] and to guarantee 
the convergence of the active contour to a useful segmentation, the gradient flow 



in (11) should be further constrained. One way to achieve this would be to use 
shape priors for T ini which would restrict the active contour to converge within 
a predefined space of expected configurations [22,24,37-42]. A conceptually new 
solution to the problem of shape priors is proposed in this paper. In particular, 
we propose to track the correct shape configuration in the manner similar to 
that used for tracking of the feature distributions. As will be shown below, the 
proposed approach is simple to implement as it requires neither registration of 
training images nor using additional computational structures apart from what 
has already been used before in the paper. An even more significant property of 
the proposed prior model is its minimally restrictive nature. The latter implies 
that the force produced by the model will be dominated by data-related forces, 
while remaining sufficiently "strong" to prevent the segmentation from diverging 
at the points of poorly observed sections of the prostate boundary. 

To construct such a "weak" prior model, we suggest to select a single param- 
eter of an expected prostate boundary T t , whose empirical pdf can be learned 
based on the same set of manually segmented images which have already been 



used to compute -Pt(z). Subsequently, the gradient flow (111 could be redefined 
in such a way that the empirical pdf of the same parameter of the active con- 
tour is forced to match, as close as possible, the learned pdf in terms of the 
Bhattacharyya coefficient ^ . 

In the current study, we use the curvature of active contours as the shape- 
tracking parameter. It is worth noting that the curvature has been chosen for 
its property of being invariant under the group of Euclidean transformations. 
Since the configuration of TRUS transducers suggests that different images of 

3 In general, an edge-detector function is defined to be strictly positive over the homoge- 
neous regions of an image, while vanishing in vicinity of the strong gradients of the image. 
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the same prostate are likely to be similar up to a rotation and/or a shift, the 
Euclidean invariance of the curvature allows one to forgo image preprocessing 
via alignment and registration. However, it should be noted that the proposed 
approach is not limited to work for the curvature alone, and thus other geometric 
descriptors of the curves could be used, should one require invariance under 
different geometric transformations [43,44]. 

Given a set of manually delineated images, the fast marching method [45] 
can be used to compute the level-set functions as the signed distance functions 
of the segmentation boundaries. Subsequently, each of the resulting level set 
functions tp(pc) can be used to compute its associated curvature according to 

«(x)=-div(^« \, (14) 



|Vp(x)| 

whose (empirical) pdf C(£ | y(x)) can in turn be defined as 



where «(x) and S e (-) are given by (|14| and (12 1, respectively. Note that the 



calculation in (14 1 provides the curvature values corresponding to all level sets 
of <^(x). Therefore, 8 £ (<p(x)) in |i"5"| ) makes the resulting pdf depend on the 
values of k(x) in immediate vicinity of the zero level set of </?(x). Finally, the 
target curvature distribution Ct(£) is obtained via averaging the curvature pdf's 
C(£ | v?( x )) over an available training images. 



The same equation (15 1 can be used to estimate the curvature pdf C(£ | 
¥>(x)) which corresponds to the level-set function </j(x) used for segmentation. 
In this case, one can measure the similarity between the target and evolving 
shapes in terms of the Bhattacharyya coefficient given as 

= / y/C t ®C(£\<p(x))d£. (16) 



Consequently, the shape priors (as defined by Ct(£)) can be incorporated into 
the segmentation process via requiring the optimal level-set function </?*( x ) to 
maximize both B(ip(x)) in ([6| and B K (ip{x)) in (161. In this case, the resulting 
optimization problem takes the form of 

p*(x) = argsup(aflfo>(x))+/3B (c ( ¥ >(x))- / <?(x) || VW(^(x))|| dxl , (17) 
v>(x) I Jn ) 

where (3 > is a shape-regularization parameter (set to be equal to 2.5 in the 
experiments reported in this paper). As a result, the optimal solution can now 
be found as a stationary point of the gradient flow defined as 

' V!X > - * e (p(x)) ( «F B (x) +div ( g( x )-^W))+pv c (x), (18) 



dr £V ^ V " \ ay ' V l|Vy>(x)| 
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where the velocity Vc(x) : f2 — > E is related to the maximization of B K (ip(x)). 
The derivation of Vc(x) is provided in Appendix B. 

Concluding this section we note that, in this paper, the gradient flow in (l8| 
was computed using the implicit discretization scheme based on the method of 
additive operator splitting (AOS) as detailed in [46]. In this case, it is common 
to perform an alternative derivation of the gradient flow, which results in an 
expression similar to (18l, with 6 e ((p(x)) replaced by ||Vy(x)||. The latter, in 
turn, could be approximated by 1, whenever y(x) is defined to be the signed 
distance function of its zero level-set [46] (which is the case in the present study) . 



3 Technical considerations 
3.1 Regularization of curvature 

As long as practical implementation of a level-set evolution is considered, it is 
standard to compute the curvature k(x) based on 

K(x) = - dlv i^(xI/ = WT^W* • (19) 

where the subscripts x and y denote partial differentiation along the corre- 
sponding directions. The partials, in turn, are normally computed by means of 
standard discretization schemes, whose numerical support is necessarily finite. 
In this case, both quantization and measurement noises can cause the discrete 
curvature to be a noisy, "jagged" function [47]. Needless to say, such behavior of 
the curvature could bias the estimation of its corresponding pdf. This situation 
can be observed in Fig. [2] Subplot A of which shows a grayscale image of the 
level-set function <£>(x) that is defined to be the signed distance function of the 
circle depicted in yellow. A zoomed fragment of the circle and of its associated 
level-set function is shown in Subplot Al, while Subplot A2 depicts the corre- 
sponding values of re(x). Note that, theoretically, the curvature of a circle is a 
constant equal to the inverse of the circle's radius. However, due to the domain 
discretization, the "discrete" circle appears to be a piecewise linear curve, and, 
as a result, the curvature of its corresponding y(x) varies considerably along 
the level-sets, as shown in Subplot A2. 

One possible way to alleviate the problem of irregularity of k(x) is to in- 



crease the numerical support of discretization of the partial derivatives in (19 1. 
Unfortunately, in the case under consideration, this approach has not resulted 
in substantial improvements. An alternative way to suppress the spurious irreg- 
ularities of discrete k(x) is to precede its computation by a regularization of the 
corresponding y(x). In particular, we suggest to regularize </?(x) via anisotropi- 
cally diffusing the latter in the direction tangent to its level-sets. Such a diffusion 
can be readily performed using the approach of [48] , which prescribes to smooth 
ip(x) through solving 

^^=div( J D(x)V^(x,r)), (20) 
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Before the regularization After the regularization 




Figure 2: (Subplot A) A level-set function and its corresponding zero level- 
set before the regularization; (Subplot Al) A zoomed fragment of the level-set 
function in Subplot A; (Subplot A2) The curvature of the level-sets in Subplot 
Al; (Subplot B) The regularized version of the level-set function in Subplot A; 
(Subplot Bl) A zoomed fragment of the regularized level-set function in Subplot 
B; (Subplot B2) The curvature of the level-sets in Subplot Bl. 
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with D £ M 2x2 being a diffusivity matrix, and r being an artificial (iteration) 
time similar to the one in ( 19 1. 

The requirement on the diffusion in ( 20 1 to propagate in the direction tangent 
to the level-sets of y(x) is controlled via a special definition of D. Specifically, 
D can be defined to have the same eigenvectors Vi and x>i as those of a smoothed 
version of the structural tensor Vip(x)Wip(x) T . Note that when the smoothing 
is negligible, the directions of v*i and V<p(x) coincide, implying v\ \\ Wip(x) 
and V2 J- Vy(x). (When the effect of smoothing cannot be ignored, the latter 
relationships will hold approximately.) Therefore, setting D = v\v^ +JV2V2, 
with 7< 1, guarantees the diffusion in (20) will smooth y(x) in the direction of 
its level-sets, while preserving the geometric shape of the zero level-set . In the 
experiments reported in this paper, the value of 7 was set to be equal to 0.1. 
Moreover, the diffusion in (20 1 was performed using the AOS-based procedure 
of [48] , with the number of iterations and the diffusion time-step being equal to 
N = 4 and At = 10, respectively. 

The value of the above procedure can be appreciated via observing Sub- 
plot B of Fig. [2] which shows regularized versions of the circular curve and its 
associated </?(x), which were previously discussed in relation to Subplot A of 
the same figure. One can see that the regularization based on (20 1 produces 
a much smoother level-set function, whose curvature better complies with its 
theoretically predicted value. 

The overall structure of the proposed method for segmentation of TRUS 
images is summarized in the pseudocode of Algorithm 1. The algorithm receives 
as its input a TRUS image u(x), the target pdf Pt(z) of a set of photometric 
features, the target pdf Ct(£) of the curvature of prostate boundary, as well 
as an initial level-set function <po(x)- At the output, the algorithm returns a 
final level-set function <pt(x-), which can be used to indicate the prostate region 
as {x I </?t(x) < 0}, or by means of its boundary {x | <^t(x) = 0}. Note that at 
Step 10 of Algorithm 1 below, the operation denoted by "AOS" refers to the 
procedure represented by equation (51) in [49]. 



3.2 Feature selection 

In the previous sections, we introduced a different way to incorporate prior 
shape information into the process of segmentation of TRUS images. It should 
be noted that, in general, the same distribution of k(x) can be shared by an 
infinite number of different shapes. For this reason, the proposed shape prior by 
itself cannot guarantee the segmentation to converge to a useful result. At the 
same time, the shape prior still contains essential information as demonstrated 
by the example of Fig. [3j The leftmost subplots of the figure show two test 
shapes, viz. "flower" and "square" , while the middle and the rightmost subplots 
show their initial and final segmentation, respectively, using the shape velocity 
Vc(x) alone. One can see that, despite the severe ambiguity of definition of the 
shape priors in the form of pdf's, the segmentations in Fig. [3] are still capable 
of converging to a close vicinity of the true shapes. 

The above examples are by no means generic, and therefore in order to 
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Test shape "Flower" Initial segmentation Segmentation after 194 iterations 




90 100 150 290 250 50 100 150 290 250 



Test shape "Square" Initial segmentation Segmentation after 104 iterations 




50 100 150 290 250 50 100 150 290 250 



Figure 3: (Upper row of subplots [left to right]) Test shape "Flower", its initial 
and final segmentations based on Vc(x) alone; (Lower row of subplots [left to 
right] ) Test shape "Square" , its initial and final segmentations based on Vc (x) 
alone. 
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Algorithm 1 Proposed segmentation procedure 

l: given: u(x), P t (z), C t (£), <A)( X ) = <p(x,* = 0) 
preset: Ar = 10, a = 0.5, /3 = 2.5, f = 
compute: {/fc(x)}^ =1 and g(x) using (13 1 (e.g. with A = 3) 
while 6 > 10~ 4 do 

Diffuse ipt (x) using ( 20 ) to result in (p t (x) 
Compute k(x) using <pt( x ) and equation ( |19[ ) 
Compute {p(z fc | yt(x))}f 3 using drb and C(£ | y t (x)) using (15) 
Compute Vb(x) using (28 1 and V"c(x) using ( [32] ) 
¥>t+i(x) <^ t (x)+At (aVs(x) + /3Vc(x)) 
¥t+i(x) «= AOS (<^ +1 (x), ff (x),Ai) 
Redistance <^ t+ i(x) by fast marching 
5 4= H^i+i - ¥>t|| 
i 4= i + 1 
end while 
return </?t(x) 



ensure the convergence of segmentation in the general case of TRUS images, 
the use of image-related features is necessary. The arsenal of possible features 
is broad, with some typical examples including (yet not limited to) Gabor and 
wavelet transform coefficients [26,50], multiresolution moments [51], local fractal 
dimension [52], and many others [15,53]. Another way to find aset of informative 
features could be to project a (large) set of arbitrary features onto a subspace 
spanned by their either independent or principal components [54]. 

Since the main contribution of this paper is the introduction of the pdf-based 
"weak" shape priors, the experimental study of the next section does not provide 
an in-depth analysis of performance of the proposed segmentation method for 
different choices of image features. The only features of TRUS images used in 
the experiments are their gray-levels as well as the values of their SRAD filtered 
versions. Despite the relative simplicity of the above choice of features, the 
segmentation results appear to be very promising as demonstrated by examples 
below. 



4 Results 

4.1 In silico experiments 

In this section, the performance of the proposed algorithm is first validated 
through an initial stage of in silico experiments, followed by subsequent in vivo 
experiments. In particular, for the case of in silico validation, the RF-images of 
prostate were simulated by convolution of synthetic reflectivity functions with a 
point spread function (PSF). The latter was obtained experimentally by imaging 
a point-target (viz. a thin steel wire in a water tank) using a single-element, 3.5 
MHz-transducer (Panametrics V383, Waltham, MA) for both transmission and 
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High contrast f4:1) Medium contrast (3:1) Low contrast {2:1) 




Figure 4: Simulated prostate images of high (4:1), medium (3:1), and low (2:1) 
contrasts. 

reception. The lateral scanning of the target was carried out mechanically with 
a lateral resolution of 0.4 mm, and the acquired RF-lincs were sampled at a rate 
of 25 MHz. The reflectivity functions, on the other hand, were generated as 2-D 
white Gaussian noise fields weighted by predefined amplitude profiles, which 
were designed to mimic the typical geometry of a prostate gland (see the first 
row of subplots in Fig. [5] for illustrative examples). Moreover, to account for 
the variability of real-life data, the profiles were subjected to random (elastic) 
deformations. 

The contrast of prostate regions was controlled by the values of the variance 
of simulated reflectivity functions inside and outside of the prostate boundaries. 
Three different ratios of the variances, viz. 4:1, 3:1, and 2:1, were used to test 
the performance of the segmentation under variable conditions. Moreover, to 
assess the robustness of the proposed method, the reflectivity functions were pre- 
multiplied by another mask to mimic the presence of shadowing artifacts as well 
as to cause the apparent boundaries of simulated prostates to be discontinuous. 
Some typical envelope images corresponding to the simulated RF-data are shown 
in Fig. |4j whose leftmost, middle, and rightmost subplots illustrate the high 
(4:1), medium (3:1), and low (2:1) contrasts, respectively. 

In the present study, 200 simulated prostate images were generated for each 
contrast category with half of them being used as training images and the rest 
used for algorithm validation. In particular, the training images together with 
their SRAD filtered versions were used to pre-compute the target pdf's -Pt(z) 
and C t (£). The test images, on the other hand, were used for segmentation, 
followed by a quantitative analysis of its performance. The latter was assessed 
in terms of the normalized mean squared error (NMSE) criterion. Specifically, 
let M be the matrix of the true binary template of a simulated prostate, and 
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Figure 5: (First row of subplots) True segmentation profiles used for simula- 
tion of prostate images; (Second row of subplots) The corresponding prostate 
images (contrast 3:1) corrupted by shadowing artifacts; (Third row of subplots) 
Segmentation results obtained by Algorithm 1 with /3 = (i.e. without shape 
priors); (Fourth row of subplots) Segmentation results obtained by Algorithm 1 
using the shape priors. 

M be its estimate. Then, the NMSE can be defined as 

where || ■ \\p denotes the Frobenius matrix norm, and E stands for the operator 
of expectation, which in the present case was approximated by a sample mean 
computed over the 100 test images. 

A subset of typical segmentation results are shown in Fig. [5] In particu- 
lar, the upper row of subplots of the figure depicts the simulated templates of 
prostate shapes used in the simulation study, while the subplots in the second 
row show the corresponding envelope images. At the first stage, these images 
were segmented by Algorithm 1 without using the "weak" shape priors (i.e. 
j3 = 0) . The results of this segmentation are demonstrated by the third row of 
subplots in Fig. [5] One can see that the discontinuity of prostate boundaries 
along with the low contrast of the shadowed regions have caused the segmen- 
tation to converge to an erroneous result. To resolve this problem, the same 
images were segmented by Algorithm 1 with the inclusion of shape priors (i.e. 
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(3 = 2.5). The resulting segmentation is shown in the bottom row of subplots of 
Fig. [5j One can see that in this case the segmentation is capable of converging 
to the true prostate shapes. 



Table 1: NMSE (± one standard deviation) of prostate segmentation 





High Contrast (4:1) 


Medium Contrast (3:1) 


Low Contrast (2:1) 


NMSE (with priors) 
NMSE (without priors) 


0.049±0.009 
0.211±0.035 


0.050±0.014 
0.247±0.048 


0.066±0.019 
0.340±0.053 



The quantitative comparison results are summarized in Table[T] which shows 
the values of NMSE obtained with and without using the shape priors for differ- 
ent contrast levels. It can be clearly seen that the inclusion of the shape priors 
results in more than 5- fold decrease in the value of NMSE. 



4.2 In vivo experiments 

In addition to the in silico experiments, in vivo experiments were conducted 
next to further test the performance of the proposed algorithm. For the latter 
case, the ultrasound images were obtained during clinical TRUS sessions using 
the Aloka 2000 ultrasound machine with a broadband 7 MHz linear transducer 
and a field of view of approximately 6 cm. Subsequently, the obtained set of 
clinical TRUS images were manually delineated by an expert radiologist. 

Apart from providing the ground truth for a comparative analysis, the man- 
ually delineated images were also used to learn the target pdf's Pt(z) and Cj(£). 
As before, the image features were defined to be the gray-level values of the 
TRUS images and their SRAD filtered versions, while the prostate morphology 
was described by the curvature of its boundary. 

In the in vivo experiments, it was found to be impossible to attain useful 
segmentation results without the use of shape priors. For this reason, only 
the prior-enhanced segmentation results are reported in this section. A typi- 
cal segmentation process is demonstrated in Fig. [6| Starting at the left and 
transitioning to the right of the figure, each column represents the results of 
initial, intermediate, and final segmentation stages. In particular, the upper 
row of subplots depicts the convergence of the active contour from its initial 
to final configuration, while Subplots B1-B3 and Subplots C1-C3 compare the 
target (red) and empirical (blue) pdf's. Specifically, the pdf's corresponding to 
the intensity values of the TRUS images (as observed within the boundaries of 
the active contour) are plotted in Subplots B1-B3. At the same time, Subplots 
C1-C3 shows the evolution of the curvature pdf's. One can see that as the seg- 
mentation converges, the appearances of the target and empirical pdf's become 
progressively more similar. 

Additional segmentation examples are demonstrated in Fig. [7] The shown 
images represent most of the principal challenges of TRUS image segmenta- 
tion, as they suffer from severe shadowing artifacts and low contrast. Moreover, 
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Figure 6: (Subplots A1-A3) Initial, intermediate, and final segmentations of a 
prostate image; (Subplots B1-B3) The target (red) and empirical (blue) pdf 's of 
the image gray-levels corresponding to the segmented regions shown in Subplots 
A1-A3, respectively; (Subplots C1-C3) The target (red) and empirical (blue) 
pdf 's of the curvature of the active contours (yellow) shown in Subplots A1-A3, 
respectively. 
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Figure 7: Examples of TRUS image segmentation: (Upper row of subplots) 
Expert manual delineations; (Lower row of subplots) Results produced by Al- 
gorithm 1. 
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the boundaries of the shown prostates are discontinuous due to the aforemen- 
tioned artifacts as well as the destructive influence of speckle noise. Despite the 
above difficulties, the proposed method yields reliable and accurate segmenta- 
tion, which closely matches the results of manual delineation. The NMSE of the 
in vivo segmentation by means of Algorithm 1 was found to be equal to 0.081 
± 0.017. 

5 Discussion and Conclusions 

A new method of automatic segmentation of TRUS images of prostate was 
proposed in this paper. The method can be categorized as a learning-based 
technique, since it takes advantage of the availability of large sets of manually 
delineated TRUS data. The latter makes it possible to learn the statistical 
properties of discrimination features characterizing the prostate region before 
the actual segmentation is carried out. Although the list of possible segmen- 
tation features is rather long, only the intensities of TRUS images and their 
SRAD-filtcred versions were used in the present study. While the resulting seg- 
mentation was found to be stable, robust, and accurate, we believe that further 
improvements are possible via a more meticulous selection of the discriminating 
features. 

The segmentation algorithm proposed in the current paper is based on the 
concept of distribution tracking, which provides an efficient and versatile frame- 
work for incorporation of an arbitrary number of statistical features in the pro- 
cess of image delineation. Moreover, one of the principal contributions of this 
research has been in introduction of a way in which the same set of ideas can be 
extended to tracking of the morphological features of prostates. This effort has 
resulted in developing the concept of "weak" priors, which are "encoded" in the 
form of pdf 's of some geometric parameters of the expected prostate boundaries. 
In this case, the segmentation is forced to converge to the shapes, which have 
the same parameters distributed similarly to the target distributions. The main 
advantage of such priors consists in their property of being minimally restric- 
tive, which suggests that the resulting estimates will have a negligible bias. On 
the other hand, the priors have been shown to be informative enough to render 
the segmentation stable and robust. 

The curvature of the prostate boundary was used as a shape descriptor to 
perform the tracking. It should be noted that, while being invariant under the 
group of Euclidean transformations, the curvature should not be used in the case 
when more complex (e.g. affine) deformations are expected. In this case, other 
geometric invariants should be considered. The definition of such invariants is 
in the focus of our current research. Among some other directions of our future 
work are an extension of the proposed method to 3-D scenarios as well as its 
verification via more extensive clinical experiments. 
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Appendix A 

To derive an expression for 5B(ip)/5tp, we first note that the first variation of 
p(zk | ¥>(x)) in d8| with respect to y(x) is equal to 



Sp(z k | ip) 
dip 



= A 



-^(xD^l^l-^-Z^x))), (22) 

where S(-) is a Dirac delta function, and 

A = I H(~tp(x))dx. (23) 
Jn 

Consequently, each of the d coefficients Bk(ip(x)) in ^ has its first variation 
defined as 

6B k (<p) 
Sip 



Pt{zk) 



Sp(z k | ip) 



2^p t {z k )p{z k | ^(x)) Sip 
%«) ^ ( S*(v(x)) - 



Pt(Zk) 



K{z k - 4(x)) dz fc 



_/ p(z fc | <p(x)) 

The above expression can be rewritten in a more concise form as 

= S (ip(x)) i (ifcMx)) - [r(z fe | p(x)) * 
095 2^ V L J Zfc =/ fc (x) 

where * stands for the operation of convolution and r(zk) is defined as 



r(z k I ip(-x)) 



Pt(zk) 



p{z k I <ys(x))' 



(24) 
(25) 

(26) 



It should be pointed out that, in practical computations, the functions r(zk \ 
<p(x)) and K(zk) are represented by their discrete values evaluated over some 
predefined subset of the real line (see Section [4] for more details) . In this case, 
the convolution in ( p5| can be efficiently computed by means of an FFT-based 
algorithm [55, Ch.8], followed by evaluating its result at the points Zk — /fe(x) 
via, e.g., either linear or cubic interpolation. 



Finally, the results in ( 25 1 can be combined together to yield the first vari- 

%>(x))Vb(x), 



ation of B(ip(x)) as given by 

SB(ip) 



Sip 



where 



1 



fe=i 



z k =I k (x) 



with 



d 

oik = Yl Bi(ip(x)). 



(27) 



(28) 



(29) 
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Appendix B 



We derive the first variation of B K (p(x)) with respect to y(x) under the as- 
sumption that ||V</?(x)|| = 1, which is valid as long as y(x) is initialized and 
sustained to be a signed-distance function, which is the case in the present study. 
Under this assumption, one can show that the first variation of C(£ | v?( x )) with 
respect to <p(x) is given by 

5C(£ hg) = A fc(y(x))IT'(g - «(x))] + ^(y(x)) - «(x)) - g(g | y(x))] 

5^3 J 5 e (ly9(x))rfx 

Consequently, the first variation of B K (<p(x)) in (Jl6| is defined as 
8B K {<p) 1 /■ / C t (0 6C(S | v») 



(30) 



Sip 



C(£ I V (x)) <V 



•d£. 



(31) 



Let L(0 = [C t (£)/C(f I ¥>(*))]' and A = /<5 £ (x)dx. Then 6B K (cp)/5(<p) be- 
comes 



jg«(y) 



1 f L( ) ^|^(x)) 

2 y mx) de - 

, ^(V(x)) 



2A 



L(0 [1T(£-k(x))-C(£Mx))K 



1 

2A 



I(f)A[i e Wx))Jf'K-K(x))],if- 



£Mx)) ( [L(0 * Jf (0] e=re(x) - B«(v(x)) 



(32) 



where * stands for the operation of convolution as before. As the next step, one 



could further expand the first (integral) term in (32). However, based on the 
results of numerical experiments, it was found that a more stable estimation of 
SB K (p)/S(p) can be achieved by directly approximating the integral using the 
trapezoidal rule together with a standard discretization of the Laplacian. Note 
that this computation is numerically efficient, as it only needs to be performed 
over the subset {x | <5 e (x) ^ 0}, which is small due to the finite support of S e (-). 
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